Analysis of Freesurfer parcellation Genetic Template, surface area and cortical thickness differences
library(easypackages)
libraries("here","heplots","ggplot2","ggseg","ggsegExtra","ggsegChen","tidyverse","patchwork","pheatmap","psych")
source(here("code","cohens_d.R"))
source(here("code","get_ggColorHue.R"))
datapath = here("data","tidy")
plotdir = here("plots")
pheno_data_gex =read.csv(file.path(datapath,"labelData_allGEXMRIsubs.csv"))
pheno_data_mri =read.csv(file.path(datapath,"labelData_all_MRI.csv"))
# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)
fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)
region_names = colnames(fs_data)[5:ncol(fs_data)]
region_names = c("CortexVol","SAtotal","CTmean",region_names)
fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")
new_fs_data = fs_data
colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
"tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
"pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use
for (icol in 1:length(region_names)){
y_var = region_names[icol]
if(is.element(y_var,c("CortexVol","SAtotal","CTmean"))){
# make formula
form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
covname2use = c("sexMale","scan_age")
}else{
# make formula
form2use = as.formula(sprintf("%s ~ Group + CortexVol + sex + scan_age", y_var))
full_model = model.matrix(~0+ Group + CortexVol + sex + scan_age, data = fs_data)
covname2use = c("CortexVol","sexMale","scan_age")
}
# run linear model
mod2use = lm(formula = form2use, data = fs_data)
# run an ANOVA
res2use = anova(mod2use)
# grab ANOVA results and store into output_res
output_res[y_var, "Fstat"] = res2use["Group", "F value"]
output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
# compute partial eta-squared as the effect size of the interaction effect
eta_sq_res = etasq(mod2use)
output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
# remove covariates before effect size computation
beta1 = mod2use$coefficients[covname2use, drop = FALSE]
beta1[is.na(beta1)] = 0
new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
# compute effect sizes
mask1 = fs_data$Group=="Good"
mask2 = fs_data$Group=="TD"
output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"]
mask1 = fs_data$Group=="Poor"
mask2 = fs_data$Group=="TD"
output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"]
mask1 = fs_data$Group=="Good"
mask2 = fs_data$Group=="Poor"
output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"]
}
output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
output_res[order(output_res$pval),]
## Fstat pval etasq
## CortexVol 18.7684243 3.563168e-08 0.0609961017
## anteromedialtemporal.lh.area 14.0080609 2.095489e-06 0.1018626300
## SAtotal 12.7918734 6.069546e-06 0.0384031239
## superiorparietal.rh.area 10.6221181 4.214024e-05 0.0406886255
## posterolateraltemporal.rh.area 10.1768804 6.297412e-05 0.0532927353
## inferiorparietal.lh.thickness 9.0672894 1.726457e-04 0.0855251058
## superiorparietal.lh.area 8.8900453 2.030244e-04 0.0469626440
## middletemporal.rh.thickness 7.2377106 9.322992e-04 0.0567858094
## dorsomedialfrontal.lh.area 6.9207431 1.252445e-03 0.0500938981
## perisylvian.lh.thickness 5.6069622 4.299179e-03 0.0280653457
## middletemporal.lh.thickness 5.1161171 6.843551e-03 0.0160809808
## inferiorparietal.rh.thickness 5.1077155 6.898359e-03 0.0452499419
## orbitofrontal.rh.area 4.7199377 9.975975e-03 0.0884379862
## dorsomedialfrontal.rh.area 4.6591244 1.057149e-02 0.0500313212
## medialprefrontal.rh.thickness 4.4135056 1.336577e-02 0.0035100067
## medialprefrontal.lh.thickness 4.3594368 1.407504e-02 0.0049069724
## inferiorparietal.rh.area 4.1241893 1.763166e-02 0.0388947624
## orbitofrontal.lh.area 3.9924278 2.000762e-02 0.0817033669
## dorsolateralprefrotal.lh.thickness 3.7760453 2.463286e-02 0.0010251563
## posterolateraltemporal.lh.area 3.6074681 2.897450e-02 0.0112156872
## occipital.lh.area 3.5466894 3.072299e-02 0.0234903844
## temporalpole.lh.thickness 3.1731609 4.407572e-02 0.0117222407
## ventromedialoccipital.lh.thickness 2.8410823 6.081898e-02 0.0428891443
## superiorparietal.rh.thickness 2.8250048 6.177622e-02 0.0223598673
## anteromedialtemporal.rh.area 2.6898314 7.045102e-02 0.0075296610
## medialtemporal.lh.thickness 2.5864636 7.790740e-02 0.0314594717
## temporalpole.rh.thickness 2.5859020 7.795001e-02 0.0019890549
## inferiorparietal.lh.area 2.5563368 8.022714e-02 0.0301583009
## CTmean 2.5204136 8.307086e-02 0.0011589273
## dorsolateralprefrotal.rh.thickness 2.4716397 8.712967e-02 0.0099715010
## ventralfrontal.rh.thickness 2.3301281 1.000286e-01 0.0204949358
## motor_premotor_SMA.rh.thickness 2.1554571 1.186455e-01 0.0549536354
## perisylvian.rh.thickness 2.1521394 1.190311e-01 0.0131780702
## superiortemporal.rh.area 2.0552479 1.308692e-01 0.0056249313
## motor_premotor_SMA.lh.thickness 1.6863034 1.879341e-01 0.0782779395
## dorsolateralprefrontal.rh.area 1.4145993 2.455442e-01 0.0080141342
## occipital.lh.thickness 1.2841529 2.792522e-01 0.0109164163
## occipital.rh.thickness 1.2733925 2.822334e-01 0.0057228488
## precuneus.rh.area 1.2492498 2.890396e-01 0.0107147043
## motorpremotor.rh.area 1.0954063 3.364843e-01 0.0149046475
## ventromedialoccipital.rh.thickness 1.0815033 3.411419e-01 0.0157302561
## parsopercularis.lh.area 1.0271582 3.599817e-01 0.0107116670
## ventralfrontal.lh.thickness 0.9478539 3.893778e-01 0.0136577437
## dorsolateralprefrontal.lh.area 0.5413956 5.828213e-01 0.0153294284
## precuneus.lh.area 0.4888912 6.140674e-01 0.0038323834
## medialtemporal.rh.thickness 0.2924523 7.467628e-01 0.0131748655
## motorpremotor.lh.area 0.2654045 7.671765e-01 0.0029461200
## occipital.rh.area 0.2082300 8.122034e-01 0.0001698168
## parsopercularis.rh.area 0.1699284 8.438520e-01 0.0001588197
## superiortemporal.lh.area 0.0788472 9.242110e-01 0.0020092765
## superiorparietal.lh.thickness 0.0213230 9.789050e-01 0.0004093216
## fdrq es.TD_vs_Good es.TD_vs_Poor
## CortexVol 1.817216e-06 -0.168185504 -0.627477606
## anteromedialtemporal.lh.area 5.343498e-05 0.073908193 -0.716130093
## SAtotal 1.031823e-04 -0.086079474 -0.471191290
## superiorparietal.rh.area 5.372881e-04 -0.239268265 0.266474769
## posterolateraltemporal.rh.area 6.423360e-04 0.166221696 -0.413272621
## inferiorparietal.lh.thickness 1.467488e-03 0.402056404 -0.376038602
## superiorparietal.lh.area 1.479178e-03 -0.338833602 0.223427618
## middletemporal.rh.thickness 5.943408e-03 -0.602325506 -0.335190443
## dorsomedialfrontal.lh.area 7.097189e-03 -0.413580399 0.117701311
## perisylvian.lh.thickness 2.192581e-02 -0.394967263 -0.059355222
## middletemporal.lh.thickness 2.931803e-02 -0.078034977 -0.363741804
## inferiorparietal.rh.thickness 2.931803e-02 0.438813968 -0.073227995
## orbitofrontal.rh.area 3.851042e-02 0.703440285 0.689021478
## dorsomedialfrontal.rh.area 3.851042e-02 -0.011856131 0.560176404
## medialprefrontal.rh.thickness 4.486419e-02 0.079704943 0.150731011
## medialprefrontal.lh.thickness 4.486419e-02 -0.005454202 0.155007839
## inferiorparietal.rh.area 5.289497e-02 -0.496123250 -0.002619246
## orbitofrontal.lh.area 5.668827e-02 0.313080777 0.746115082
## dorsolateralprefrotal.lh.thickness 6.611979e-02 0.056890925 0.081060057
## posterolateraltemporal.lh.area 7.388497e-02 0.074128125 -0.187712338
## occipital.lh.area 7.461298e-02 0.395346960 0.163865789
## temporalpole.lh.thickness 1.021755e-01 -0.011829570 -0.257701161
## ventromedialoccipital.lh.thickness 1.312745e-01 -0.415214100 -0.555129522
## superiorparietal.rh.thickness 1.312745e-01 0.308387021 -0.015645109
## anteromedialtemporal.rh.area 1.437201e-01 -0.067557401 -0.230844671
## medialtemporal.lh.thickness 1.460901e-01 0.100526613 0.450751114
## temporalpole.rh.thickness 1.460901e-01 -0.102905181 -0.020598602
## inferiorparietal.lh.area 1.460901e-01 -0.297574970 -0.448702546
## CTmean 1.460901e-01 0.022296408 0.086644973
## dorsolateralprefrotal.rh.thickness 1.481204e-01 0.140067779 -0.103150861
## ventralfrontal.rh.thickness 1.645632e-01 -0.323190343 -0.295450988
## motor_premotor_SMA.rh.thickness 1.839571e-01 0.313202751 0.608448832
## perisylvian.rh.thickness 1.839571e-01 -0.119473742 -0.292404471
## superiortemporal.rh.area 1.963038e-01 0.182570796 0.138964153
## motor_premotor_SMA.lh.thickness 2.738468e-01 0.573830868 0.746481899
## dorsolateralprefrontal.rh.area 3.478543e-01 0.167982081 -0.040570128
## occipital.lh.thickness 3.779749e-01 -0.111820856 0.148339120
## occipital.rh.thickness 3.779749e-01 -0.013600907 0.166091886
## precuneus.rh.area 3.779749e-01 -0.117983093 0.134026057
## motorpremotor.rh.area 4.243472e-01 -0.207271672 -0.303815461
## ventromedialoccipital.rh.thickness 4.243472e-01 -0.302852519 -0.232545925
## parsopercularis.lh.area 4.371207e-01 0.026319646 -0.249586445
## ventralfrontal.lh.thickness 4.618202e-01 -0.264409577 -0.259488718
## dorsolateralprefrontal.lh.area 6.755429e-01 0.171882416 0.316840463
## precuneus.lh.area 6.959431e-01 -0.077260554 0.076319613
## medialtemporal.rh.thickness 8.279326e-01 0.089540153 0.324278245
## motorpremotor.lh.area 8.324682e-01 -0.010298772 0.112012146
## occipital.rh.area 8.629661e-01 0.030998595 0.028420544
## parsopercularis.rh.area 8.782949e-01 -0.027358771 -0.026814437
## superiortemporal.lh.area 9.426953e-01 -0.114992084 -0.083632069
## superiorparietal.lh.thickness 9.789050e-01 0.023147126 0.057755463
## es.Good_vs_Poor tstat.TD_vs_Good
## CortexVol -0.416871265 -0.92763624
## anteromedialtemporal.lh.area -0.758745394 0.37243309
## SAtotal -0.372984599 -0.53069497
## superiorparietal.rh.area 0.547055816 -1.56200620
## posterolateraltemporal.rh.area -0.604517309 1.15445234
## inferiorparietal.lh.thickness -0.744128402 2.52598615
## superiorparietal.lh.area 0.528719718 -2.00216978
## middletemporal.rh.thickness 0.292282835 -3.22118075
## dorsomedialfrontal.lh.area 0.538454727 -2.36785475
## perisylvian.lh.thickness 0.344432195 -2.31179712
## middletemporal.lh.thickness -0.226876400 -0.59981901
## inferiorparietal.rh.thickness -0.523112641 2.46021689
## orbitofrontal.rh.area -0.022666089 3.81615464
## dorsomedialfrontal.rh.area 0.481823664 -0.03266780
## medialprefrontal.rh.thickness 0.075467838 0.62367574
## medialprefrontal.lh.thickness 0.148087974 -0.10676434
## inferiorparietal.rh.area 0.396157491 -2.97475639
## orbitofrontal.lh.area 0.472178056 1.54212952
## dorsolateralprefrotal.lh.thickness 0.021731886 0.08491989
## posterolateraltemporal.lh.area -0.267730933 0.36341567
## occipital.lh.area -0.219531034 2.06315786
## temporalpole.lh.thickness -0.229334676 -0.05302691
## ventromedialoccipital.lh.thickness -0.054958989 -2.36766848
## superiorparietal.rh.thickness -0.350452936 1.98361558
## anteromedialtemporal.rh.area -0.144985360 -0.48454790
## medialtemporal.lh.thickness 0.348994398 0.52586747
## temporalpole.rh.thickness 0.082064561 -0.68884090
## inferiorparietal.lh.area -0.177660163 -1.59491429
## CTmean 0.059512817 0.06023806
## dorsolateralprefrotal.rh.thickness -0.269073157 0.28318744
## ventralfrontal.rh.thickness 0.049662712 -2.02068983
## motor_premotor_SMA.rh.thickness 0.317960692 1.78735887
## perisylvian.rh.thickness -0.180757516 -0.92800055
## superiortemporal.rh.area -0.041254182 1.26918573
## motor_premotor_SMA.lh.thickness 0.132796743 3.26084740
## dorsolateralprefrontal.rh.area -0.225924020 1.12326898
## occipital.lh.thickness 0.273712261 -0.76184410
## occipital.rh.thickness 0.172930032 0.02467738
## precuneus.rh.area 0.288530088 -0.69967304
## motorpremotor.rh.area -0.101723418 -1.08822596
## ventromedialoccipital.rh.thickness 0.095756438 -1.51414221
## parsopercularis.lh.area -0.235320467 0.06273803
## ventralfrontal.lh.thickness -0.008185191 -1.58638447
## dorsolateralprefrontal.lh.area 0.153830025 1.12998643
## precuneus.lh.area 0.161333323 -0.29115018
## medialtemporal.rh.thickness 0.193392174 0.62787481
## motorpremotor.lh.area 0.122913072 -0.03900672
## occipital.rh.area 0.001588241 -0.08532334
## parsopercularis.rh.area 0.003390924 0.01719381
## superiortemporal.lh.area 0.030124858 -0.30438030
## superiorparietal.lh.thickness 0.026537911 0.53721444
## tstat.TD_vs_Poor tstat.Good_vs_Poor
## CortexVol -3.059262385 2.51015442
## anteromedialtemporal.lh.area -3.515199051 4.06483461
## SAtotal -2.131564929 2.19003479
## superiorparietal.rh.area 1.984471136 -2.77322638
## posterolateraltemporal.rh.area -2.047477385 3.54666956
## inferiorparietal.lh.thickness -2.131633301 4.08238560
## superiorparietal.lh.area 1.551859981 -2.78194967
## middletemporal.rh.thickness -1.778064339 -1.60535282
## dorsomedialfrontal.lh.area 0.509988831 -3.41232980
## perisylvian.lh.thickness -0.296281168 -2.01527339
## middletemporal.lh.thickness -1.040308667 1.32359190
## inferiorparietal.rh.thickness -0.867470744 2.38090177
## orbitofrontal.rh.area 3.758467841 0.06798017
## dorsomedialfrontal.rh.area 3.017287570 -2.63515835
## medialprefrontal.rh.thickness 0.871125758 -0.31661320
## medialprefrontal.lh.thickness 0.819719437 -0.78027216
## inferiorparietal.rh.area 0.499151365 -1.80659810
## orbitofrontal.lh.area 3.924457099 -3.01607583
## dorsolateralprefrotal.lh.thickness 0.722697259 -0.12019995
## posterolateraltemporal.lh.area -0.612513896 1.73587140
## occipital.lh.area 0.228185483 0.88449841
## temporalpole.lh.thickness -1.353811398 1.51173501
## ventromedialoccipital.lh.thickness -2.908615860 -0.02696405
## superiorparietal.rh.thickness -0.149573963 1.97211361
## anteromedialtemporal.rh.area -0.893081728 0.90645129
## medialtemporal.lh.thickness 2.275557563 -1.70721752
## temporalpole.rh.thickness -0.065063096 -0.25861361
## inferiorparietal.lh.area -2.228747825 1.12256572
## CTmean 0.003809356 -0.36930388
## dorsolateralprefrotal.rh.thickness -0.083335625 1.44315454
## ventralfrontal.rh.thickness -1.470130934 -0.01116362
## motor_premotor_SMA.rh.thickness 3.062774577 -1.85472165
## perisylvian.rh.thickness -1.210572588 0.71872792
## superiortemporal.rh.area 0.356212566 0.34908797
## motor_premotor_SMA.lh.thickness 3.933261727 -0.91827562
## dorsolateralprefrontal.rh.area -0.331176983 1.28323113
## occipital.lh.thickness 0.315153460 -1.87591783
## occipital.rh.thickness 0.828635994 -0.78044635
## precuneus.rh.area 1.160264079 -1.74582006
## motorpremotor.rh.area -1.779218478 0.37142742
## ventromedialoccipital.rh.thickness -1.134943163 -0.54892205
## parsopercularis.lh.area -1.564937351 1.22474572
## ventralfrontal.lh.thickness -1.209652541 0.50358676
## dorsolateralprefrontal.lh.area 1.193483568 -0.85504869
## precuneus.lh.area 0.844134840 -0.70718029
## medialtemporal.rh.thickness 1.638318531 -0.72499110
## motorpremotor.lh.area 0.837840504 -0.34973418
## occipital.rh.area -0.333265662 -0.45488366
## parsopercularis.rh.area -0.860511040 -0.10119179
## superiortemporal.lh.area -0.762767692 -0.15225809
## superiorparietal.lh.thickness -0.143367881 -0.20504914
## pval.TD_vs_Good pval.TD_vs_Poor
## CortexVol 0.3553559051 0.0026877095
## anteromedialtemporal.lh.area 0.7101961158 0.0006040535
## SAtotal 0.5965569661 0.0348955186
## superiorparietal.rh.area 0.1207947235 0.0492929914
## posterolateraltemporal.rh.area 0.2504996384 0.0426100630
## inferiorparietal.lh.thickness 0.0127767912 0.0349040119
## superiorparietal.lh.area 0.0474136371 0.1231083165
## middletemporal.rh.thickness 0.0016246125 0.0777132171
## dorsomedialfrontal.lh.area 0.0194129637 0.6109174538
## perisylvian.lh.thickness 0.0224120292 0.7674845535
## middletemporal.lh.thickness 0.5497040259 0.3001127370
## inferiorparietal.rh.thickness 0.0152403869 0.3872703845
## orbitofrontal.rh.area 0.0002113866 0.0002564376
## dorsomedialfrontal.rh.area 0.9739911935 0.0030651648
## medialprefrontal.rh.thickness 0.5339679967 0.3852791198
## medialprefrontal.lh.thickness 0.9151457399 0.4138649853
## inferiorparietal.rh.area 0.0035153952 0.6185099704
## orbitofrontal.lh.area 0.1255499164 0.0001396919
## dorsolateralprefrotal.lh.thickness 0.9324598646 0.4711541563
## posterolateraltemporal.lh.area 0.7169034100 0.5412595960
## occipital.lh.area 0.0411514717 0.8198578416
## temporalpole.lh.thickness 0.9577944449 0.1781276166
## ventromedialoccipital.lh.thickness 0.0194223111 0.0042658904
## superiorparietal.rh.thickness 0.0494727774 0.8813307850
## anteromedialtemporal.rh.area 0.6288383274 0.3734510978
## medialtemporal.lh.thickness 0.5999046095 0.0244961519
## temporalpole.rh.thickness 0.4921898421 0.9482229469
## inferiorparietal.lh.area 0.1132367114 0.0275351297
## CTmean 0.9520607209 0.9969663319
## dorsolateralprefrotal.rh.thickness 0.7774978971 0.9337118827
## ventralfrontal.rh.thickness 0.0454314672 0.1439240155
## motor_premotor_SMA.rh.thickness 0.0762835306 0.0026621984
## perisylvian.rh.thickness 0.3551815765 0.2282388836
## superiortemporal.rh.area 0.2067142418 0.7222544192
## motor_premotor_SMA.lh.thickness 0.0014289517 0.0001351941
## dorsolateralprefrontal.rh.area 0.2634589483 0.7410399734
## occipital.lh.thickness 0.4475767309 0.7531466001
## occipital.rh.thickness 0.9803513384 0.4088176208
## precuneus.rh.area 0.4854212466 0.2480514687
## motorpremotor.rh.area 0.2785725200 0.0775229375
## ventromedialoccipital.rh.thickness 0.1324942477 0.2584714928
## parsopercularis.lh.area 0.9500745288 0.1200101061
## ventralfrontal.lh.thickness 0.1151586014 0.2285906765
## dorsolateralprefrontal.lh.area 0.2606284515 0.2348369502
## precuneus.lh.area 0.7714151924 0.4001329423
## medialtemporal.rh.thickness 0.5312222760 0.1037549440
## motorpremotor.lh.area 0.9689467714 0.4036463233
## occipital.rh.area 0.9321397658 0.7394665581
## parsopercularis.rh.area 0.9863091943 0.3910795681
## superiortemporal.lh.area 0.7613403104 0.4469732857
## superiorparietal.lh.thickness 0.5920669333 0.8862199452
## pval.Good_vs_Poor fdrq.TD_vs_Good
## CortexVol 1.334511e-02 0.69704428
## anteromedialtemporal.lh.area 8.475100e-05 0.96215984
## SAtotal 3.037760e-02 0.87414672
## superiorparietal.rh.area 6.408005e-03 0.35564245
## posterolateraltemporal.rh.area 5.514425e-04 0.58419158
## inferiorparietal.lh.thickness 7.930211e-05 0.12381723
## superiorparietal.lh.area 6.248267e-03 0.19408551
## middletemporal.rh.thickness 1.109600e-01 0.02761841
## dorsomedialfrontal.lh.area 8.703711e-04 0.12381723
## perisylvian.lh.thickness 4.603818e-02 0.12700150
## middletemporal.lh.thickness 1.880737e-01 0.87414672
## inferiorparietal.rh.thickness 1.879256e-02 0.12381723
## orbitofrontal.rh.area 9.459109e-01 0.01078072
## dorsomedialfrontal.rh.area 9.482980e-03 0.98630919
## medialprefrontal.rh.thickness 7.520695e-01 0.87414672
## medialprefrontal.lh.thickness 4.367176e-01 0.98630919
## inferiorparietal.rh.area 7.325004e-02 0.04482129
## orbitofrontal.lh.area 3.107251e-03 0.35564245
## dorsolateralprefrotal.lh.thickness 9.045193e-01 0.98630919
## posterolateraltemporal.lh.area 8.507058e-02 0.96215984
## occipital.lh.area 3.781391e-01 0.19408551
## temporalpole.lh.thickness 1.331457e-01 0.98630919
## ventromedialoccipital.lh.thickness 9.785318e-01 0.12381723
## superiorparietal.rh.thickness 5.082244e-02 0.19408551
## anteromedialtemporal.rh.area 3.664551e-01 0.89085430
## medialtemporal.lh.thickness 9.028439e-02 0.87414672
## temporalpole.rh.thickness 7.963622e-01 0.86557524
## inferiorparietal.lh.area 2.637909e-01 0.35564245
## CTmean 7.125261e-01 0.98630919
## dorsolateralprefrotal.rh.thickness 1.514984e-01 0.96713153
## ventralfrontal.rh.thickness 9.911108e-01 0.19408551
## motor_premotor_SMA.rh.thickness 6.601136e-02 0.27789000
## perisylvian.rh.thickness 4.736608e-01 0.69704428
## superiortemporal.rh.area 7.276154e-01 0.52712132
## motor_premotor_SMA.lh.thickness 3.602573e-01 0.02761841
## dorsolateralprefrontal.rh.area 2.018042e-01 0.58419158
## occipital.lh.thickness 6.301777e-02 0.84542271
## occipital.rh.thickness 4.366155e-01 0.98630919
## precuneus.rh.area 8.331901e-02 0.86557524
## motorpremotor.rh.area 7.109531e-01 0.59196660
## ventromedialoccipital.rh.thickness 5.840460e-01 0.35564245
## parsopercularis.lh.area 2.229926e-01 0.98630919
## ventralfrontal.lh.thickness 6.154453e-01 0.35564245
## dorsolateralprefrontal.lh.area 3.941733e-01 0.58419158
## precuneus.lh.area 4.807814e-01 0.96713153
## medialtemporal.rh.thickness 4.698234e-01 0.87414672
## motorpremotor.lh.area 7.271315e-01 0.98630919
## occipital.rh.area 6.499883e-01 0.98630919
## parsopercularis.rh.area 9.195617e-01 0.98630919
## superiortemporal.lh.area 8.792308e-01 0.96713153
## superiorparietal.lh.thickness 8.378699e-01 0.87414672
## fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## CortexVol 0.022331915 0.075622317
## anteromedialtemporal.lh.area 0.007701682 0.002161151
## SAtotal 0.148342050 0.140841621
## superiorparietal.rh.area 0.179567326 0.046686893
## posterolateraltemporal.rh.area 0.167162555 0.009374523
## inferiorparietal.lh.thickness 0.148342050 0.002161151
## superiorparietal.lh.area 0.330448639 0.046686893
## middletemporal.rh.thickness 0.247710879 0.282947979
## dorsomedialfrontal.lh.area 0.788600212 0.011097232
## perisylvian.lh.thickness 0.869815827 0.195662284
## middletemporal.lh.thickness 0.566879614 0.417033003
## inferiorparietal.rh.thickness 0.603060407 0.095842051
## orbitofrontal.rh.area 0.004359440 0.984519528
## dorsomedialfrontal.rh.area 0.022331915 0.060453999
## medialprefrontal.rh.thickness 0.603060407 0.891989414
## medialprefrontal.lh.thickness 0.603060407 0.696018673
## inferiorparietal.rh.area 0.788600212 0.233484491
## orbitofrontal.lh.area 0.003562144 0.031693961
## dorsolateralprefrotal.lh.thickness 0.649428702 0.977034258
## posterolateraltemporal.lh.area 0.726427353 0.241033306
## occipital.lh.area 0.908972824 0.665003166
## temporalpole.lh.thickness 0.432595640 0.323353722
## ventromedialoccipital.lh.thickness 0.027195051 0.991110843
## superiorparietal.rh.thickness 0.941608692 0.199380340
## anteromedialtemporal.rh.area 0.603060407 0.665003166
## medialtemporal.lh.thickness 0.138811527 0.242342319
## temporalpole.rh.thickness 0.967187406 0.923056146
## inferiorparietal.lh.area 0.140429161 0.517435932
## CTmean 0.996966332 0.883532933
## dorsolateralprefrotal.rh.thickness 0.967187406 0.351200902
## ventralfrontal.rh.thickness 0.367006240 0.991110843
## motor_premotor_SMA.rh.thickness 0.022331915 0.224438630
## perisylvian.rh.thickness 0.499028519 0.700567122
## superiortemporal.rh.area 0.869815827 0.883532933
## motor_premotor_SMA.lh.thickness 0.003562144 0.665003166
## dorsolateralprefrontal.rh.area 0.869815827 0.428833975
## occipital.lh.thickness 0.869815827 0.224438630
## occipital.rh.thickness 0.603060407 0.696018673
## precuneus.rh.area 0.506024996 0.241033306
## motorpremotor.rh.area 0.247710879 0.883532933
## ventromedialoccipital.rh.thickness 0.507001774 0.827398532
## parsopercularis.lh.area 0.330448639 0.454904910
## ventralfrontal.lh.thickness 0.499028519 0.848316434
## dorsolateralprefrontal.lh.area 0.499028519 0.670094525
## precuneus.lh.area 0.603060407 0.700567122
## medialtemporal.rh.thickness 0.311264832 0.700567122
## motorpremotor.lh.area 0.603060407 0.883532933
## occipital.rh.area 0.869815827 0.872352757
## parsopercularis.rh.area 0.603060407 0.977034258
## superiortemporal.lh.area 0.633212155 0.974799407
## superiorparietal.lh.thickness 0.941608692 0.949585942
regions2use = colnames(fs_data)[5:(ncol(fs_data)-4)]
full_regnames = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
"middle temporal","perisylvian","dorsolateral prefrontal","occipital",
"temporal pole","superior parietal","medial prefrontal","medial temporal",
"motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
"middle temporal","perisylvian","dorsolateral prefrontal","occipital",
"temporal pole","superior parietal","medial prefrontal","medial temporal")
sa_ordering = c(1,12,7,10,
4,6,9,3,
8,11,2,5)
ct_ordering = c(1,7,6,3,
10,4,11,5,
8,2,12,9)
gensim_order = c(sa_ordering,sa_ordering, ct_ordering,ct_ordering)
df4plot = output_res[regions2use,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")]
df4plot$phenotype = "CT"
df4plot$phenotype[1:24] = "SA"
df4plot$phenotype = factor(df4plot$phenotype)
df4plot$hemisphere = "RH"
df4plot$hemisphere[1:12] = "LH"
df4plot$hemisphere[25:36] = "LH"
df4plot$hemisphere = factor(df4plot$hemisphere)
df4plot$full_region_names = factor(full_regnames)
df4plot$gensim_gradient = gensim_order
ct_df = subset(df4plot, df4plot$phenotype=="CT")
sa_df = subset(df4plot, df4plot$phenotype=="SA")
ct_df$cluster = "dorsal"
ct_df$cluster[c(6,7,8,9,10,11,12)] = "ventral"
ct_df$cluster = factor(ct_df$cluster)
sa_df$cluster = "anterior"
sa_df$cluster[c(6,7,8,9,10,11,12)] = "posterior"
sa_df$cluster = factor(sa_df$cluster)
dotSize = 10
fontSize = 20
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Good)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.TD_vs_Good
## t = -1.838, df = 22, p-value = 0.07961
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.66968163 0.04520194
## sample estimates:
## cor
## -0.3648475
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Good by cluster
## t = 0.55381, df = 20.769, p-value = 0.5856
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1553724 0.2680582
## sample estimates:
## mean in group dorsal mean in group ventral
## 0.007270414 -0.049072461
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.TD_vs_Poor
## t = -0.88864, df = 22, p-value = 0.3838
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5483668 0.2348869
## sample estimates:
## cor
## -0.1861478
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Poor by cluster
## t = -1.1479, df = 18.319, p-value = 0.2658
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3867892 0.1132490
## sample estimates:
## mean in group dorsal mean in group ventral
## -0.05450477 0.08226532
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.Good_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.Good_vs_Poor
## t = 0.80097, df = 22, p-value = 0.4317
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2521901 0.5353739
## sample estimates:
## cor
## 0.1683315
p = ggplot(data = ct_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.Good_vs_Poor by cluster
## t = -1.7772, df = 15.521, p-value = 0.09513
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.41512987 0.03702278
## sample estimates:
## mean in group dorsal mean in group ventral
## -0.05560184 0.13345171
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Good)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.TD_vs_Good
## t = -0.84267, df = 22, p-value = 0.4085
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.541591 0.243973
## sample estimates:
## cor
## -0.1768267
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Good by cluster
## t = 1.1527, df = 15.531, p-value = 0.2665
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1017404 0.3429604
## sample estimates:
## mean in group anterior mean in group posterior
## 0.03042461 -0.09018539
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.TD_vs_Poor
## t = -0.82792, df = 22, p-value = 0.4166
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5393989 0.2468832
## sample estimates:
## cor
## -0.173825
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Poor by cluster
## t = 1.6011, df = 10.357, p-value = 0.1394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09637807 0.59688583
## sample estimates:
## mean in group anterior mean in group posterior
## 0.1092268 -0.1410271
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.Good_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.Good_vs_Poor
## t = -0.13004, df = 22, p-value = 0.8977
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4263442 0.3799303
## sample estimates:
## cor
## -0.02771394
p = ggplot(data = sa_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.Good_vs_Poor by cluster
## t = 0.64854, df = 9.4989, p-value = 0.532
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2769209 0.5020282
## sample estimates:
## mean in group anterior mean in group posterior
## 0.07141322 -0.04114044
dotSize = 6
fontSize=28
fdr_thresh = 0.01
output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
output_res_fdr[order(output_res_fdr$pval),]
## Fstat pval etasq fdrq
## CortexVol 18.768424 3.563168e-08 0.06099610 1.817216e-06
## anteromedialtemporal.lh.area 14.008061 2.095489e-06 0.10186263 5.343498e-05
## SAtotal 12.791873 6.069546e-06 0.03840312 1.031823e-04
## superiorparietal.rh.area 10.622118 4.214024e-05 0.04068863 5.372881e-04
## posterolateraltemporal.rh.area 10.176880 6.297412e-05 0.05329274 6.423360e-04
## inferiorparietal.lh.thickness 9.067289 1.726457e-04 0.08552511 1.467488e-03
## superiorparietal.lh.area 8.890045 2.030244e-04 0.04696264 1.479178e-03
## middletemporal.rh.thickness 7.237711 9.322992e-04 0.05678581 5.943408e-03
## dorsomedialfrontal.lh.area 6.920743 1.252445e-03 0.05009390 7.097189e-03
## es.TD_vs_Good es.TD_vs_Poor es.Good_vs_Poor
## CortexVol -0.16818550 -0.6274776 -0.4168713
## anteromedialtemporal.lh.area 0.07390819 -0.7161301 -0.7587454
## SAtotal -0.08607947 -0.4711913 -0.3729846
## superiorparietal.rh.area -0.23926826 0.2664748 0.5470558
## posterolateraltemporal.rh.area 0.16622170 -0.4132726 -0.6045173
## inferiorparietal.lh.thickness 0.40205640 -0.3760386 -0.7441284
## superiorparietal.lh.area -0.33883360 0.2234276 0.5287197
## middletemporal.rh.thickness -0.60232551 -0.3351904 0.2922828
## dorsomedialfrontal.lh.area -0.41358040 0.1177013 0.5384547
## tstat.TD_vs_Good tstat.TD_vs_Poor
## CortexVol -0.9276362 -3.0592624
## anteromedialtemporal.lh.area 0.3724331 -3.5151991
## SAtotal -0.5306950 -2.1315649
## superiorparietal.rh.area -1.5620062 1.9844711
## posterolateraltemporal.rh.area 1.1544523 -2.0474774
## inferiorparietal.lh.thickness 2.5259861 -2.1316333
## superiorparietal.lh.area -2.0021698 1.5518600
## middletemporal.rh.thickness -3.2211807 -1.7780643
## dorsomedialfrontal.lh.area -2.3678548 0.5099888
## tstat.Good_vs_Poor pval.TD_vs_Good
## CortexVol 2.510154 0.355355905
## anteromedialtemporal.lh.area 4.064835 0.710196116
## SAtotal 2.190035 0.596556966
## superiorparietal.rh.area -2.773226 0.120794723
## posterolateraltemporal.rh.area 3.546670 0.250499638
## inferiorparietal.lh.thickness 4.082386 0.012776791
## superiorparietal.lh.area -2.781950 0.047413637
## middletemporal.rh.thickness -1.605353 0.001624612
## dorsomedialfrontal.lh.area -3.412330 0.019412964
## pval.TD_vs_Poor pval.Good_vs_Poor
## CortexVol 0.0026877095 1.334511e-02
## anteromedialtemporal.lh.area 0.0006040535 8.475100e-05
## SAtotal 0.0348955186 3.037760e-02
## superiorparietal.rh.area 0.0492929914 6.408005e-03
## posterolateraltemporal.rh.area 0.0426100630 5.514425e-04
## inferiorparietal.lh.thickness 0.0349040119 7.930211e-05
## superiorparietal.lh.area 0.1231083165 6.248267e-03
## middletemporal.rh.thickness 0.0777132171 1.109600e-01
## dorsomedialfrontal.lh.area 0.6109174538 8.703711e-04
## fdrq.TD_vs_Good fdrq.TD_vs_Poor
## CortexVol 0.69704428 0.022331915
## anteromedialtemporal.lh.area 0.96215984 0.007701682
## SAtotal 0.87414672 0.148342050
## superiorparietal.rh.area 0.35564245 0.179567326
## posterolateraltemporal.rh.area 0.58419158 0.167162555
## inferiorparietal.lh.thickness 0.12381723 0.148342050
## superiorparietal.lh.area 0.19408551 0.330448639
## middletemporal.rh.thickness 0.02761841 0.247710879
## dorsomedialfrontal.lh.area 0.12381723 0.788600212
## fdrq.Good_vs_Poor
## CortexVol 0.075622317
## anteromedialtemporal.lh.area 0.002161151
## SAtotal 0.140841621
## superiorparietal.rh.area 0.046686893
## posterolateraltemporal.rh.area 0.009374523
## inferiorparietal.lh.thickness 0.002161151
## superiorparietal.lh.area 0.046686893
## middletemporal.rh.thickness 0.282947979
## dorsomedialfrontal.lh.area 0.011097232
df2plot = output_res[4:nrow(output_res),]
df2plot$region = rownames(df2plot)
new_region_labels = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor-supplementary motor area","ventral frontal",
"ventromedial occipital","inferior parietal","middle temporal","perisylvian",
"dorsolateral prefrontal","occipital","temporal pole","superior parietal",
"medial prefrontal","medial temporal",
"motor-premotor-supplementary motor area","ventral frontal",
"ventromedial occipital","inferior parietal","middle temporal","perisylvian",
"dorsolateral prefrontal","occipital","temporal pole","superior parietal",
"medial prefrontal","medial temporal")
hemi_labels = c("LH","LH","LH","LH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"RH","RH","RH","RH")
metric_labels = c("SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT")
df2plot$region = new_region_labels
df2plot$hemisphere = hemi_labels
df2plot$metric = metric_labels
colors2use = get_ggColorHue(7)
# SA F-stat
d2plot = subset(df2plot, df2plot$metric=="SA" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))
# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot,
display_numbers = round(mat2plot,digits=2),
number_color = "black", fontsize_number = 18,
show_rownames=TRUE,
labels_col = colnames(mat2plot),
color = colorRampPalette(c('light blue','white','red'))(100),
cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
c("#785e00","#c49a00","#ffcc12"),
c("#306800","#53b400","#76ff01"),
c("#007459","#00c094","#0effc8"),
c("#007b9f","#00b6eb","#39d2ff"),
c("#6a3eff","#a58aff","#e0d7ff"),
c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(0.6,1.3), c(0.8,1.15), c(0.75,1.35), c(0.8,1.35),
c(0.8,1.4), c(0.8, 1.15), c(0.7,1.1))
for (ireg in 1:length(sig_reg)){
region = sig_reg[ireg]
reg_name = reg_names[ireg]
metric_name = metric_names[ireg]
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
p1 = p1 + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
print(p1)
tmp_df = d2plot[region,]
if (tmp_df$hemisphere=="LH"){
HEMI = "left"
} else if (tmp_df$hemisphere=="RH"){
HEMI = "right"
}
p = ggseg(.data = tmp_df,
atlas = chenAr,
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere=HEMI) +
scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) +
theme(text = element_text(size=fontSize))
print(p)
}










cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CortexVol"
reg_name = "Total Cortical Volume"
metric_name = "Volume"
ylims = c(450000,700000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p1

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "SAtotal"
reg_name = "Total Surface Area"
metric_name = "SA"
ylims = c(90000,220000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p1

# SA es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# SA es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# SA es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT F-stat
d2plot = subset(df2plot, df2plot$metric=="CT" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))
# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot,
display_numbers = round(mat2plot,digits=2),
number_color = "black", fontsize_number = 18,
show_rownames=TRUE,
labels_col = colnames(mat2plot),
color = colorRampPalette(c('light blue','white','red'))(100),
cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_manual(values = colors2use[1:4]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_manual(values = colors2use[5:7]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
c("#785e00","#c49a00","#ffcc12"),
c("#306800","#53b400","#76ff01"),
c("#007459","#00c094","#0effc8"),
c("#007b9f","#00b6eb","#39d2ff"),
c("#6a3eff","#a58aff","#e0d7ff"),
c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(-0.4,0.1), c(-0.4,0.4), c(-0.4,0.3), c(0,0.7),
c(-0.3,0.3), c(-0.2,0.5), c(0,0.8))
for (ireg in 1:length(sig_reg)){
region = sig_reg[ireg]
reg_name = reg_names[ireg]
metric_name = metric_names[ireg]
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
p1 = p1 + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
print(p1)
tmp_df = d2plot[region,]
if (tmp_df$hemisphere=="LH"){
HEMI = "left"
} else if (tmp_df$hemisphere=="RH"){
HEMI = "right"
}
p = ggseg(.data = tmp_df,
atlas = chenTh,
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere=HEMI) +
scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) +
theme(text = element_text(size=fontSize))
print(p)
}




# CT es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

Do modeling of local differences without adjustment for total cortical volume.
# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)
fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)
region_names = colnames(fs_data)[5:ncol(fs_data)]
region_names = c("CortexVol","SAtotal","CTmean",region_names)
fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")
new_fs_data = fs_data
colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
"tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
"pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use
for (icol in 1:length(region_names)){
y_var = region_names[icol]
# make formula
form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
covname2use = c("sexMale","scan_age")
# run linear model
mod2use = lm(formula = form2use, data = fs_data)
# run an ANOVA
res2use = anova(mod2use)
# grab ANOVA results and store into output_res
output_res[y_var, "Fstat"] = res2use["Group", "F value"]
output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
# compute partial eta-squared as the effect size of the interaction effect
eta_sq_res = etasq(mod2use)
output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
# remove covariates before effect size computation
beta1 = mod2use$coefficients[covname2use, drop = FALSE]
beta1[is.na(beta1)] = 0
new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
# compute effect sizes
mask1 = fs_data$Group=="Good"
mask2 = fs_data$Group=="TD"
output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"]
mask1 = fs_data$Group=="Poor"
mask2 = fs_data$Group=="TD"
output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"]
mask1 = fs_data$Group=="Good"
mask2 = fs_data$Group=="Poor"
output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])
d2use = subset(fs_data, mask1 | mask2)
mod2use = lm(formula = form2use, data = d2use)
res = summary(mod2use)
output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"]
output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"]
}
output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
output_res[order(output_res$pval),]
## Fstat pval etasq
## CortexVol 18.76842425 3.563168e-08 6.099610e-02
## anteromedialtemporal.lh.area 13.96198552 2.171932e-06 1.187908e-01
## SAtotal 12.79187344 6.069546e-06 3.840312e-02
## superiorparietal.rh.area 9.97848940 7.518206e-05 6.112751e-02
## posterolateraltemporal.rh.area 9.79945558 8.843829e-05 7.466123e-02
## inferiorparietal.lh.thickness 9.11412690 1.650865e-04 8.726838e-02
## superiorparietal.lh.area 8.39947078 3.178795e-04 6.424037e-02
## middletemporal.rh.thickness 7.25613913 9.152796e-04 5.855042e-02
## dorsomedialfrontal.lh.area 6.95235733 1.214613e-03 4.948907e-02
## perisylvian.lh.thickness 5.52500358 4.641913e-03 2.703457e-02
## inferiorparietal.rh.thickness 5.13340268 6.727702e-03 4.630736e-02
## middletemporal.lh.thickness 4.88755620 8.498974e-03 3.322598e-02
## orbitofrontal.rh.area 4.73445787 9.833290e-03 8.651774e-02
## dorsomedialfrontal.rh.area 4.67781875 1.037894e-02 4.928688e-02
## medialprefrontal.rh.thickness 4.26035165 1.546793e-02 1.189294e-02
## medialprefrontal.lh.thickness 4.12360606 1.763387e-02 1.682454e-02
## inferiorparietal.rh.area 4.05447261 1.884306e-02 4.356279e-02
## orbitofrontal.lh.area 3.97768109 2.028472e-02 7.394201e-02
## occipital.lh.area 3.54678869 3.071012e-02 2.442052e-02
## posterolateraltemporal.lh.area 3.47096883 3.304108e-02 2.333083e-02
## dorsolateralprefrotal.lh.thickness 3.38761197 3.581075e-02 1.222666e-02
## temporalpole.lh.thickness 3.10692393 4.698360e-02 2.235669e-02
## ventromedialoccipital.lh.thickness 2.84207875 6.074741e-02 4.976358e-02
## superiorparietal.rh.thickness 2.80424833 6.302173e-02 2.572353e-02
## anteromedialtemporal.rh.area 2.70365990 6.949654e-02 8.384856e-03
## medialtemporal.lh.thickness 2.57867720 7.848677e-02 2.562891e-02
## CTmean 2.52041363 8.307086e-02 1.158927e-03
## inferiorparietal.lh.area 2.51632609 8.340244e-02 2.073572e-02
## dorsolateralprefrotal.rh.thickness 2.43567482 9.022542e-02 6.603310e-03
## temporalpole.rh.thickness 2.25922675 1.071864e-01 9.199317e-03
## perisylvian.rh.thickness 2.16026648 1.180743e-01 1.186715e-02
## motor_premotor_SMA.rh.thickness 2.10997052 1.240297e-01 4.047085e-02
## superiortemporal.rh.area 2.06577098 1.295134e-01 5.537640e-03
## ventralfrontal.rh.thickness 1.98772991 1.398011e-01 1.206705e-02
## motor_premotor_SMA.lh.thickness 1.61382458 2.017999e-01 5.746284e-02
## dorsolateralprefrontal.rh.area 1.38860475 2.519036e-01 1.261047e-02
## occipital.lh.thickness 1.28537276 2.789040e-01 1.422125e-02
## occipital.rh.thickness 1.26440340 2.847364e-01 1.084307e-02
## precuneus.rh.area 1.25553820 2.872390e-01 1.159653e-02
## motorpremotor.rh.area 1.09409666 3.369095e-01 2.026216e-02
## ventromedialoccipital.rh.thickness 1.08302969 3.406168e-01 1.771466e-02
## parsopercularis.lh.area 1.02483692 3.607994e-01 7.390424e-03
## ventralfrontal.lh.thickness 0.83696514 4.345882e-01 6.276432e-03
## dorsolateralprefrontal.lh.area 0.53335542 5.874954e-01 8.638171e-03
## precuneus.lh.area 0.48186041 6.183743e-01 8.298034e-03
## medialtemporal.rh.thickness 0.27931431 7.566073e-01 3.548265e-03
## motorpremotor.lh.area 0.26497128 7.675066e-01 1.317870e-03
## occipital.rh.area 0.20422404 8.154556e-01 2.499252e-03
## parsopercularis.rh.area 0.17081341 8.431061e-01 1.638071e-04
## superiortemporal.lh.area 0.07917468 9.239085e-01 2.428987e-03
## superiorparietal.lh.thickness 0.02132311 9.789049e-01 2.247971e-05
## fdrq es.TD_vs_Good es.TD_vs_Poor
## CortexVol 1.817216e-06 -0.168185504 -0.627477606
## anteromedialtemporal.lh.area 5.538426e-05 0.058564872 -0.779840965
## SAtotal 1.031823e-04 -0.086079474 -0.471191290
## superiorparietal.rh.area 9.020706e-04 -0.192365758 0.422061004
## posterolateraltemporal.rh.area 9.020706e-04 0.132202091 -0.530611464
## inferiorparietal.lh.thickness 1.403235e-03 0.403204612 -0.371601402
## superiorparietal.lh.area 2.315979e-03 -0.284174733 0.382993616
## middletemporal.rh.thickness 5.834907e-03 -0.607526883 -0.368109847
## dorsomedialfrontal.lh.area 6.882808e-03 -0.417777838 0.101394689
## perisylvian.lh.thickness 2.367376e-02 -0.417755510 -0.141804387
## inferiorparietal.rh.thickness 3.119207e-02 0.436965414 -0.080961775
## middletemporal.lh.thickness 3.612064e-02 -0.111703511 -0.484911933
## orbitofrontal.rh.area 3.780900e-02 0.693881222 0.663368868
## dorsomedialfrontal.rh.area 3.780900e-02 -0.017052068 0.536607213
## medialprefrontal.rh.thickness 5.259096e-02 0.108471803 0.272612686
## medialprefrontal.lh.thickness 5.620796e-02 0.034413715 0.303600784
## inferiorparietal.rh.area 5.652917e-02 -0.465146831 0.088250958
## orbitofrontal.lh.area 5.747337e-02 0.302307914 0.684921768
## occipital.lh.area 8.243242e-02 0.405773203 0.206264627
## posterolateraltemporal.lh.area 8.425477e-02 0.040680967 -0.313873302
## dorsolateralprefrotal.lh.thickness 8.696896e-02 0.100397423 0.287923700
## temporalpole.lh.thickness 1.089165e-01 -0.036104455 -0.356907842
## ventromedialoccipital.lh.thickness 1.339212e-01 -0.426319170 -0.598130176
## superiorparietal.rh.thickness 1.339212e-01 0.287207325 -0.083702106
## anteromedialtemporal.rh.area 1.417729e-01 -0.068797641 -0.236087202
## medialtemporal.lh.thickness 1.519116e-01 0.086741514 0.393647416
## CTmean 1.519116e-01 0.022296408 0.086644973
## inferiorparietal.lh.area 1.519116e-01 -0.264343885 -0.361967769
## dorsolateralprefrotal.rh.thickness 1.586723e-01 0.158394449 -0.019007650
## temporalpole.rh.thickness 1.822168e-01 -0.155643965 -0.252523004
## perisylvian.rh.thickness 1.942512e-01 -0.112868711 -0.270616154
## motor_premotor_SMA.rh.thickness 1.976724e-01 0.283654397 0.508152496
## superiortemporal.rh.area 2.001571e-01 0.181205381 0.133269647
## ventralfrontal.rh.thickness 2.097016e-01 -0.245560738 -0.034176376
## motor_premotor_SMA.lh.thickness 2.940513e-01 0.532800679 0.594555794
## dorsolateralprefrontal.rh.area 3.568634e-01 0.141028172 -0.131031873
## occipital.lh.thickness 3.756202e-01 -0.101511168 0.187461609
## occipital.rh.thickness 3.756202e-01 0.003233587 0.234784901
## precuneus.rh.area 3.756202e-01 -0.116016165 0.142181444
## motorpremotor.rh.area 4.236940e-01 -0.219240608 -0.350517881
## ventromedialoccipital.rh.thickness 4.236940e-01 -0.309481377 -0.275857079
## parsopercularis.lh.area 4.381135e-01 0.038454375 -0.187192775
## ventralfrontal.lh.thickness 5.154418e-01 -0.201740570 -0.035264557
## dorsolateralprefrontal.lh.area 6.809605e-01 0.146391233 0.230406377
## precuneus.lh.area 7.008242e-01 -0.055471505 0.161218164
## medialtemporal.rh.thickness 8.328263e-01 0.054942375 0.168683420
## motorpremotor.lh.area 8.328263e-01 -0.023422791 0.062490351
## occipital.rh.area 8.664216e-01 0.056561428 0.118308443
## parsopercularis.rh.area 8.775186e-01 -0.027403626 -0.027013701
## superiortemporal.lh.area 9.423867e-01 -0.120637855 -0.104382282
## superiorparietal.lh.thickness 9.789049e-01 0.011660023 0.008979405
## es.Good_vs_Poor tstat.TD_vs_Good
## CortexVol -0.416871265 -0.92763624
## anteromedialtemporal.lh.area -0.794803988 0.16733131
## SAtotal -0.372984599 -0.53069497
## superiorparietal.rh.area 0.627981093 -1.28721394
## posterolateraltemporal.rh.area -0.684868006 0.96665266
## inferiorparietal.lh.thickness -0.741115244 2.47252195
## superiorparietal.lh.area 0.617690683 -1.58921178
## middletemporal.rh.thickness 0.267713645 -3.33001733
## dorsomedialfrontal.lh.area 0.525846171 -2.35106489
## perisylvian.lh.thickness 0.272324607 -2.37407644
## inferiorparietal.rh.thickness -0.529004839 2.51701196
## middletemporal.lh.thickness -0.320912774 -0.64723145
## orbitofrontal.rh.area -0.043209574 3.72905434
## dorsomedialfrontal.rh.area 0.466219000 -0.04161523
## medialprefrontal.rh.thickness 0.161735686 0.88959051
## medialprefrontal.lh.thickness 0.247253521 0.14898278
## inferiorparietal.rh.area 0.446921552 -2.83757132
## orbitofrontal.lh.area 0.418212137 1.63638457
## occipital.lh.area -0.188622586 2.13738216
## posterolateraltemporal.lh.area -0.360688505 0.09481595
## dorsolateralprefrotal.lh.thickness 0.173019231 0.46757198
## temporalpole.lh.thickness -0.300282103 -0.24817897
## ventromedialoccipital.lh.thickness -0.085302255 -2.38119144
## superiorparietal.rh.thickness -0.403949755 1.77629274
## anteromedialtemporal.rh.area -0.148614009 -0.53739563
## medialtemporal.lh.thickness 0.303841318 0.45893648
## CTmean 0.059512817 0.06023806
## inferiorparietal.lh.area -0.116200804 -1.32874444
## dorsolateralprefrotal.rh.thickness -0.196139656 0.47428576
## temporalpole.rh.thickness -0.080405337 -0.92866917
## perisylvian.rh.thickness -0.163250435 -0.71061831
## motor_premotor_SMA.rh.thickness 0.238945563 1.65793388
## superiortemporal.rh.area -0.045635530 1.32733836
## ventralfrontal.rh.thickness 0.218413802 -1.59784005
## motor_premotor_SMA.lh.thickness 0.034411163 3.10024557
## dorsolateralprefrontal.rh.area -0.297885547 0.89205477
## occipital.lh.thickness 0.305012767 -0.66781162
## occipital.rh.thickness 0.219594316 0.02112775
## precuneus.rh.area 0.295072659 -0.64623092
## motorpremotor.rh.area -0.139743512 -1.13212189
## ventromedialoccipital.rh.thickness 0.068451005 -1.66827098
## parsopercularis.lh.area -0.195904940 0.12501934
## ventralfrontal.lh.thickness 0.142797772 -1.35421352
## dorsolateralprefrontal.lh.area 0.088096781 0.92201234
## precuneus.lh.area 0.222805118 -0.19224401
## medialtemporal.rh.thickness 0.086489168 0.38406302
## motorpremotor.lh.area 0.086177916 -0.24000071
## occipital.rh.area 0.071541806 0.13971882
## parsopercularis.rh.area 0.003238283 -0.09907955
## superiortemporal.lh.area 0.016663201 -0.21163657
## superiorparietal.lh.thickness -0.003378940 0.27808249
## tstat.TD_vs_Poor tstat.Good_vs_Poor
## CortexVol -3.059262385 2.51015442
## anteromedialtemporal.lh.area -3.457718886 4.44548917
## SAtotal -2.131564929 2.19003479
## superiorparietal.rh.area 2.459248096 -3.58518048
## posterolateraltemporal.rh.area -2.687567084 4.10526246
## inferiorparietal.lh.thickness -2.045168291 4.13088855
## superiorparietal.lh.area 2.027310188 -3.41805335
## middletemporal.rh.thickness -1.747812746 -1.47491113
## dorsomedialfrontal.lh.area 0.493902033 -3.20521339
## perisylvian.lh.thickness -0.804089909 -1.53042696
## inferiorparietal.rh.thickness -1.006601493 2.65923159
## middletemporal.lh.thickness -2.249001377 1.86010176
## orbitofrontal.rh.area 4.159885029 0.32847869
## dorsomedialfrontal.rh.area 3.102459847 -2.54332578
## medialprefrontal.rh.thickness 1.200505533 -0.82509176
## medialprefrontal.lh.thickness 1.480004889 -1.40025248
## inferiorparietal.rh.area 0.661789105 -2.43272506
## orbitofrontal.lh.area 3.660090310 -2.40201886
## occipital.lh.area 0.791523855 0.99772299
## posterolateraltemporal.lh.area -1.116709014 2.19083651
## dorsolateralprefrotal.lh.thickness 1.511150424 -0.92199454
## temporalpole.lh.thickness -1.870933210 1.82402736
## ventromedialoccipital.lh.thickness -3.457098374 0.21658321
## superiorparietal.rh.thickness -0.406187547 2.20945034
## anteromedialtemporal.rh.area -0.952204003 0.83690241
## medialtemporal.lh.thickness 2.192865606 -1.44986930
## CTmean 0.003809356 -0.36930388
## inferiorparietal.lh.area -2.167118278 0.72881076
## dorsolateralprefrotal.rh.thickness 0.019698750 0.99601950
## temporalpole.rh.thickness -1.230129724 0.75667370
## perisylvian.rh.thickness -1.517818628 0.77574387
## motor_premotor_SMA.rh.thickness 2.832075777 -1.28000273
## superiortemporal.rh.area 0.112341163 0.34644699
## ventralfrontal.rh.thickness -0.549063137 -1.31651197
## motor_premotor_SMA.lh.thickness 3.517773482 -0.12326995
## dorsolateralprefrontal.rh.area -0.610137269 1.68378954
## occipital.lh.thickness 0.546421993 -1.93570525
## occipital.rh.thickness 1.235174967 -1.27724201
## precuneus.rh.area 1.048263005 -1.84337987
## motorpremotor.rh.area -2.128265711 0.67938285
## ventromedialoccipital.rh.thickness -0.804993424 -0.33280224
## parsopercularis.lh.area -1.396669225 0.99693457
## ventralfrontal.lh.thickness -0.236705025 -0.77527381
## dorsolateralprefrontal.lh.area 0.943125104 -0.51787823
## precuneus.lh.area 1.148868563 -1.27347487
## medialtemporal.rh.thickness 1.351096705 -0.04991803
## motorpremotor.lh.area 0.669118895 -0.41733396
## occipital.rh.area 0.260126462 -0.55310296
## parsopercularis.rh.area -0.598475549 -0.13375612
## superiortemporal.lh.area -1.371031861 -0.10594187
## superiorparietal.lh.thickness 0.180689829 -0.04619990
## pval.TD_vs_Good pval.TD_vs_Poor
## CortexVol 0.3553559051 0.0026877095
## anteromedialtemporal.lh.area 0.8673755167 0.0007336348
## SAtotal 0.5965569661 0.0348955186
## superiorparietal.rh.area 0.2003607435 0.0152169867
## posterolateraltemporal.rh.area 0.3355548864 0.0081243665
## inferiorparietal.lh.thickness 0.0147388815 0.0428253247
## superiorparietal.lh.area 0.1144990012 0.0446439975
## middletemporal.rh.thickness 0.0011370894 0.0828217214
## dorsomedialfrontal.lh.area 0.0202592759 0.6221962228
## perisylvian.lh.thickness 0.0190909801 0.4227907979
## inferiorparietal.rh.thickness 0.0130803984 0.3159679753
## middletemporal.lh.thickness 0.5186503403 0.0261692628
## orbitofrontal.rh.area 0.0002883625 0.0000569994
## dorsomedialfrontal.rh.area 0.9668707826 0.0023477810
## medialprefrontal.rh.thickness 0.3753677892 0.2320930412
## medialprefrontal.lh.thickness 0.8818035363 0.1412547090
## inferiorparietal.rh.area 0.0052934622 0.5092600271
## orbitofrontal.lh.area 0.1042349448 0.0003634346
## occipital.lh.area 0.0344834500 0.4300583771
## posterolateraltemporal.lh.area 0.9246104541 0.2661473395
## dorsolateralprefrotal.lh.thickness 0.6408920214 0.1331404116
## temporalpole.lh.thickness 0.8043967759 0.0635682953
## ventromedialoccipital.lh.thickness 0.0187420175 0.0007351844
## superiorparietal.rh.thickness 0.0780791851 0.6852625412
## anteromedialtemporal.rh.area 0.5919347149 0.3427338726
## medialtemporal.lh.thickness 0.6470642744 0.0300690708
## CTmean 0.9520607209 0.9969663319
## inferiorparietal.lh.area 0.1863147768 0.0320206760
## dorsolateralprefrotal.rh.thickness 0.6361105984 0.9843134327
## temporalpole.rh.thickness 0.3548222583 0.2208363467
## perisylvian.rh.thickness 0.4786237488 0.1314515201
## motor_premotor_SMA.rh.thickness 0.0997988784 0.0053501067
## superiortemporal.rh.area 0.1867779598 0.9107235653
## ventralfrontal.rh.thickness 0.1125637166 0.5838897264
## motor_premotor_SMA.lh.thickness 0.0023821153 0.0005973643
## dorsolateralprefrontal.rh.area 0.3740506397 0.5428198559
## occipital.lh.thickness 0.5054656127 0.5856983221
## occipital.rh.thickness 0.9831769052 0.2189593711
## precuneus.rh.area 0.5192958825 0.2964330144
## motorpremotor.rh.area 0.2597162243 0.0351732624
## ventromedialoccipital.rh.thickness 0.0977256551 0.4222710686
## parsopercularis.lh.area 0.9007059164 0.1648569587
## ventralfrontal.lh.thickness 0.1780729303 0.8132525724
## dorsolateralprefrontal.lh.area 0.3582703860 0.3473398905
## precuneus.lh.area 0.8478581339 0.2526877319
## medialtemporal.rh.thickness 0.7015740123 0.1789757279
## motorpremotor.lh.area 0.8107166551 0.5045882392
## occipital.rh.area 0.8891034877 0.7951715449
## parsopercularis.rh.area 0.9212313305 0.5505482448
## superiortemporal.lh.area 0.8327298285 0.1726915381
## superiorparietal.lh.thickness 0.7814012556 0.8568881688
## pval.Good_vs_Poor fdrq.TD_vs_Good
## CortexVol 1.334511e-02 0.70902805
## anteromedialtemporal.lh.area 1.911862e-05 0.98239861
## SAtotal 3.037760e-02 0.91667439
## superiorparietal.rh.area 4.814179e-04 0.51091990
## posterolateraltemporal.rh.area 7.238701e-05 0.70902805
## inferiorparietal.lh.thickness 6.564083e-05 0.11480256
## superiorparietal.lh.area 8.518760e-04 0.36496557
## middletemporal.rh.thickness 1.427507e-01 0.02899578
## dorsomedialfrontal.lh.area 1.713282e-03 0.11480256
## perisylvian.lh.thickness 1.284379e-01 0.11480256
## inferiorparietal.rh.thickness 8.857225e-03 0.11480256
## middletemporal.lh.thickness 6.522152e-02 0.85432548
## orbitofrontal.rh.area 7.430994e-01 0.01470649
## dorsomedialfrontal.rh.area 1.219960e-02 0.98317691
## medialprefrontal.rh.thickness 4.108913e-01 0.70902805
## medialprefrontal.lh.thickness 1.639150e-01 0.98239861
## inferiorparietal.rh.area 1.640132e-02 0.06749164
## orbitofrontal.lh.area 1.777626e-02 0.36496557
## occipital.lh.area 3.203412e-01 0.17586560
## posterolateraltemporal.lh.area 3.031810e-02 0.98239861
## dorsolateralprefrotal.lh.thickness 3.583076e-01 0.91667439
## temporalpole.lh.thickness 7.053658e-02 0.98239861
## ventromedialoccipital.lh.thickness 8.288862e-01 0.11480256
## superiorparietal.rh.thickness 2.896481e-02 0.36200349
## anteromedialtemporal.rh.area 4.042446e-01 0.91667439
## medialtemporal.lh.thickness 1.495994e-01 0.91667439
## CTmean 7.125261e-01 0.98317691
## inferiorparietal.lh.area 4.674807e-01 0.50135137
## dorsolateralprefrotal.rh.thickness 3.211649e-01 0.91667439
## temporalpole.rh.thickness 4.506692e-01 0.70902805
## perisylvian.rh.thickness 4.393654e-01 0.85432548
## motor_premotor_SMA.rh.thickness 2.029146e-01 0.36496557
## superiortemporal.rh.area 7.295893e-01 0.50135137
## ventralfrontal.rh.thickness 1.904111e-01 0.36496557
## motor_premotor_SMA.lh.thickness 9.020913e-01 0.04049596
## dorsolateralprefrontal.rh.area 9.471804e-02 0.70902805
## occipital.lh.thickness 5.516167e-02 0.85432548
## occipital.rh.thickness 2.038841e-01 0.98317691
## precuneus.rh.area 6.764219e-02 0.85432548
## motorpremotor.rh.area 4.981514e-01 0.63073940
## ventromedialoccipital.rh.thickness 7.398411e-01 0.36496557
## parsopercularis.lh.area 3.207222e-01 0.98239861
## ventralfrontal.lh.thickness 4.396421e-01 0.50135137
## dorsolateralprefrontal.lh.area 6.054586e-01 0.70902805
## precuneus.lh.area 2.052126e-01 0.98239861
## medialtemporal.rh.thickness 9.602674e-01 0.96703445
## motorpremotor.lh.area 6.771501e-01 0.98239861
## occipital.rh.area 5.811808e-01 0.98239861
## parsopercularis.rh.area 8.938105e-01 0.98239861
## superiortemporal.lh.area 9.157983e-01 0.98239861
## superiorparietal.lh.thickness 9.632247e-01 0.98239861
## fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## CortexVol 0.019581884 0.0756223169
## anteromedialtemporal.lh.area 0.007498881 0.0009750497
## SAtotal 0.119589092 0.1106612739
## superiorparietal.rh.area 0.077606632 0.0061380785
## posterolateraltemporal.rh.area 0.046038077 0.0012305792
## inferiorparietal.lh.thickness 0.133931993 0.0012305792
## superiorparietal.lh.area 0.133931993 0.0086891350
## middletemporal.rh.thickness 0.222310936 0.3466803416
## dorsomedialfrontal.lh.area 0.721181985 0.0145628942
## perisylvian.lh.thickness 0.592783168 0.3275166530
## inferiorparietal.rh.thickness 0.503573961 0.0645312131
## middletemporal.lh.thickness 0.119589092 0.1998536554
## orbitofrontal.rh.area 0.002906969 0.8421793049
## dorsomedialfrontal.rh.area 0.019581884 0.0756223169
## medialprefrontal.rh.thickness 0.422740896 0.6350138435
## medialprefrontal.lh.thickness 0.327454098 0.3634636481
## inferiorparietal.rh.area 0.665955420 0.0824171823
## orbitofrontal.lh.area 0.007498881 0.0824171823
## occipital.lh.area 0.592783168 0.5459802676
## posterolateraltemporal.lh.area 0.452450477 0.1106612739
## dorsolateralprefrotal.lh.thickness 0.323341000 0.5894737178
## temporalpole.lh.thickness 0.180110170 0.1998536554
## ventromedialoccipital.lh.thickness 0.007498881 0.9189825748
## superiorparietal.rh.thickness 0.776630880 0.1106612739
## anteromedialtemporal.rh.area 0.521009836 0.6350138435
## medialtemporal.lh.thickness 0.119589092 0.3467986259
## CTmean 0.996966332 0.8421793049
## inferiorparietal.lh.area 0.119589092 0.6443652932
## dorsolateralprefrotal.rh.thickness 0.996966332 0.5459802676
## temporalpole.rh.thickness 0.417135322 0.6384480199
## perisylvian.rh.thickness 0.323341000 0.6384480199
## motor_premotor_SMA.rh.thickness 0.034106930 0.3876238194
## superiortemporal.rh.area 0.947895956 0.8421793049
## ventralfrontal.rh.thickness 0.694665452 0.3876238194
## motor_premotor_SMA.lh.thickness 0.007498881 0.9531778237
## dorsolateralprefrontal.rh.area 0.684828305 0.2542431500
## occipital.lh.thickness 0.694665452 0.1875496803
## occipital.rh.thickness 0.417135322 0.3876238194
## precuneus.rh.area 0.487680120 0.1998536554
## motorpremotor.rh.area 0.119589092 0.6685716490
## ventromedialoccipital.rh.thickness 0.592783168 0.8421793049
## parsopercularis.lh.area 0.365110485 0.5459802676
## ventralfrontal.lh.thickness 0.882465557 0.6384480199
## dorsolateralprefrontal.lh.area 0.521009836 0.7719597195
## precuneus.lh.area 0.444381873 0.3876238194
## medialtemporal.rh.thickness 0.365110485 0.9632246542
## motorpremotor.lh.area 0.665955420 0.8421793049
## occipital.rh.area 0.881603235 0.7600056174
## parsopercularis.rh.area 0.684828305 0.9531778237
## superiortemporal.lh.area 0.365110485 0.9531778237
## superiorparietal.lh.thickness 0.910443679 0.9632246542
regions2use = colnames(fs_data)[5:(ncol(fs_data)-4)]
full_regnames = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
"middle temporal","perisylvian","dorsolateral prefrontal","occipital",
"temporal pole","superior parietal","medial prefrontal","medial temporal",
"motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
"middle temporal","perisylvian","dorsolateral prefrontal","occipital",
"temporal pole","superior parietal","medial prefrontal","medial temporal")
sa_ordering = c(1,12,7,10,
4,6,9,3,
8,11,2,5)
ct_ordering = c(1,7,6,3,
10,4,11,5,
8,2,12,9)
gensim_order = c(sa_ordering,sa_ordering, ct_ordering,ct_ordering)
df4plot = output_res[regions2use,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")]
df4plot$phenotype = "CT"
df4plot$phenotype[1:24] = "SA"
df4plot$phenotype = factor(df4plot$phenotype)
df4plot$hemisphere = "RH"
df4plot$hemisphere[1:12] = "LH"
df4plot$hemisphere[25:36] = "LH"
df4plot$hemisphere = factor(df4plot$hemisphere)
df4plot$full_region_names = factor(full_regnames)
df4plot$gensim_gradient = gensim_order
ct_df = subset(df4plot, df4plot$phenotype=="CT")
sa_df = subset(df4plot, df4plot$phenotype=="SA")
ct_df$cluster = "dorsal"
ct_df$cluster[c(6,7,8,9,10,11,12)] = "ventral"
ct_df$cluster = factor(ct_df$cluster)
sa_df$cluster = "anterior"
sa_df$cluster[c(6,7,8,9,10,11,12)] = "posterior"
sa_df$cluster = factor(sa_df$cluster)
dotSize = 10
fontSize = 20
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Good)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.TD_vs_Good
## t = -1.661, df = 22, p-value = 0.1109
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.64971696 0.08041565
## sample estimates:
## cor
## -0.3338096
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Good ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Good by cluster
## t = 0.50429, df = 19.576, p-value = 0.6197
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1627292 0.2663074
## sample estimates:
## mean in group dorsal mean in group ventral
## 0.005766454 -0.046022637
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.TD_vs_Poor
## t = -0.26326, df = 22, p-value = 0.7948
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4492788 0.3553922
## sample estimates:
## cor
## -0.05603835
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Poor ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Poor by cluster
## t = -1.2347, df = 13.777, p-value = 0.2376
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4439168 0.1198526
## sample estimates:
## mean in group dorsal mean in group ventral
## -0.06447487 0.09755724
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("CT ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.Good_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: ct_df$gensim_gradient and ct_df$es.Good_vs_Poor
## t = 1.3098, df = 22, p-value = 0.2038
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1507977 0.6065479
## sample estimates:
## cor
## 0.2689602
p = ggplot(data = ct_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Dorsal-Ventral Clusters ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.Good_vs_Poor ~ cluster, data = ct_df)
##
## Welch Two Sample t-test
##
## data: es.Good_vs_Poor by cluster
## t = -1.8651, df = 14.85, p-value = 0.08205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.44196835 0.02964677
## sample estimates:
## mean in group dorsal mean in group ventral
## -0.0636193 0.1425415
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Good)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.TD_vs_Good
## t = -0.5984, df = 22, p-value = 0.5557
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5042089 0.2917377
## sample estimates:
## cor
## -0.1265525
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters TD vs ASD Good") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Good ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Good by cluster
## t = 1.197, df = 15.567, p-value = 0.2492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09268998 0.33188479
## sample estimates:
## mean in group anterior mean in group posterior
## 0.03176589 -0.08783151
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.TD_vs_Poor
## t = -0.13074, df = 22, p-value = 0.8972
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4264660 0.3798029
## sample estimates:
## cor
## -0.02786268
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters TD vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Poor ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.TD_vs_Poor by cluster
## t = 1.5331, df = 11.265, p-value = 0.1529
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1049842 0.5914930
## sample estimates:
## mean in group anterior mean in group posterior
## 0.1089166 -0.1343378
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
geom_point(size=dotSize, colour="black", pch=21) +
geom_smooth(method=lm, colour="black") + xlim(0,12) +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
guides(colour=FALSE, fill=FALSE) +
xlab("Brain Region") +
ylab("Effect Size (Cohen's d)") +
ggtitle("SA ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.Good_vs_Poor)
##
## Pearson's product-moment correlation
##
## data: sa_df$gensim_gradient and sa_df$es.Good_vs_Poor
## t = 0.32019, df = 22, p-value = 0.7518
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3447631 0.4588948
## sample estimates:
## cor
## 0.0681055
p = ggplot(data = sa_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
# scale_fill_manual(values = c("#1e90ff","#ff8d1e")) +
scale_fill_manual(values = c("blue","red")) +
guides(colour=FALSE, fill=FALSE) +
xlab("Clusters") +
ylab("Effect Size (Cohen's d)") +
ggtitle("Anterior-Posterior Clusters ASD Good vs ASD Poor") +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.Good_vs_Poor ~ cluster, data = sa_df)
##
## Welch Two Sample t-test
##
## data: es.Good_vs_Poor by cluster
## t = 0.5723, df = 10.21, p-value = 0.5795
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2949289 0.4995674
## sample estimates:
## mean in group anterior mean in group posterior
## 0.06610517 -0.03621407
dotSize = 6
fontSize=28
fdr_thresh = 0.01
output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
output_res_fdr[order(output_res_fdr$pval),]
## Fstat pval etasq fdrq
## CortexVol 18.768424 3.563168e-08 0.06099610 1.817216e-06
## anteromedialtemporal.lh.area 13.961986 2.171932e-06 0.11879079 5.538426e-05
## SAtotal 12.791873 6.069546e-06 0.03840312 1.031823e-04
## superiorparietal.rh.area 9.978489 7.518206e-05 0.06112751 9.020706e-04
## posterolateraltemporal.rh.area 9.799456 8.843829e-05 0.07466123 9.020706e-04
## inferiorparietal.lh.thickness 9.114127 1.650865e-04 0.08726838 1.403235e-03
## superiorparietal.lh.area 8.399471 3.178795e-04 0.06424037 2.315979e-03
## middletemporal.rh.thickness 7.256139 9.152796e-04 0.05855042 5.834907e-03
## dorsomedialfrontal.lh.area 6.952357 1.214613e-03 0.04948907 6.882808e-03
## es.TD_vs_Good es.TD_vs_Poor es.Good_vs_Poor
## CortexVol -0.16818550 -0.6274776 -0.4168713
## anteromedialtemporal.lh.area 0.05856487 -0.7798410 -0.7948040
## SAtotal -0.08607947 -0.4711913 -0.3729846
## superiorparietal.rh.area -0.19236576 0.4220610 0.6279811
## posterolateraltemporal.rh.area 0.13220209 -0.5306115 -0.6848680
## inferiorparietal.lh.thickness 0.40320461 -0.3716014 -0.7411152
## superiorparietal.lh.area -0.28417473 0.3829936 0.6176907
## middletemporal.rh.thickness -0.60752688 -0.3681098 0.2677136
## dorsomedialfrontal.lh.area -0.41777784 0.1013947 0.5258462
## tstat.TD_vs_Good tstat.TD_vs_Poor
## CortexVol -0.9276362 -3.059262
## anteromedialtemporal.lh.area 0.1673313 -3.457719
## SAtotal -0.5306950 -2.131565
## superiorparietal.rh.area -1.2872139 2.459248
## posterolateraltemporal.rh.area 0.9666527 -2.687567
## inferiorparietal.lh.thickness 2.4725220 -2.045168
## superiorparietal.lh.area -1.5892118 2.027310
## middletemporal.rh.thickness -3.3300173 -1.747813
## dorsomedialfrontal.lh.area -2.3510649 0.493902
## tstat.Good_vs_Poor pval.TD_vs_Good
## CortexVol 2.510154 0.355355905
## anteromedialtemporal.lh.area 4.445489 0.867375517
## SAtotal 2.190035 0.596556966
## superiorparietal.rh.area -3.585180 0.200360743
## posterolateraltemporal.rh.area 4.105262 0.335554886
## inferiorparietal.lh.thickness 4.130889 0.014738882
## superiorparietal.lh.area -3.418053 0.114499001
## middletemporal.rh.thickness -1.474911 0.001137089
## dorsomedialfrontal.lh.area -3.205213 0.020259276
## pval.TD_vs_Poor pval.Good_vs_Poor
## CortexVol 0.0026877095 1.334511e-02
## anteromedialtemporal.lh.area 0.0007336348 1.911862e-05
## SAtotal 0.0348955186 3.037760e-02
## superiorparietal.rh.area 0.0152169867 4.814179e-04
## posterolateraltemporal.rh.area 0.0081243665 7.238701e-05
## inferiorparietal.lh.thickness 0.0428253247 6.564083e-05
## superiorparietal.lh.area 0.0446439975 8.518760e-04
## middletemporal.rh.thickness 0.0828217214 1.427507e-01
## dorsomedialfrontal.lh.area 0.6221962228 1.713282e-03
## fdrq.TD_vs_Good fdrq.TD_vs_Poor
## CortexVol 0.70902805 0.019581884
## anteromedialtemporal.lh.area 0.98239861 0.007498881
## SAtotal 0.91667439 0.119589092
## superiorparietal.rh.area 0.51091990 0.077606632
## posterolateraltemporal.rh.area 0.70902805 0.046038077
## inferiorparietal.lh.thickness 0.11480256 0.133931993
## superiorparietal.lh.area 0.36496557 0.133931993
## middletemporal.rh.thickness 0.02899578 0.222310936
## dorsomedialfrontal.lh.area 0.11480256 0.721181985
## fdrq.Good_vs_Poor
## CortexVol 0.0756223169
## anteromedialtemporal.lh.area 0.0009750497
## SAtotal 0.1106612739
## superiorparietal.rh.area 0.0061380785
## posterolateraltemporal.rh.area 0.0012305792
## inferiorparietal.lh.thickness 0.0012305792
## superiorparietal.lh.area 0.0086891350
## middletemporal.rh.thickness 0.3466803416
## dorsomedialfrontal.lh.area 0.0145628942
df2plot = output_res[4:nrow(output_res),]
df2plot$region = rownames(df2plot)
new_region_labels = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor","occipital","posterolateral temporal","superior parietal",
"orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
"anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
"motor-premotor-supplementary motor area","ventral frontal",
"ventromedial occipital","inferior parietal","middle temporal","perisylvian",
"dorsolateral prefrontal","occipital","temporal pole","superior parietal",
"medial prefrontal","medial temporal",
"motor-premotor-supplementary motor area","ventral frontal",
"ventromedial occipital","inferior parietal","middle temporal","perisylvian",
"dorsolateral prefrontal","occipital","temporal pole","superior parietal",
"medial prefrontal","medial temporal")
hemi_labels = c("LH","LH","LH","LH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"LH","LH","LH","LH",
"RH","RH","RH","RH",
"RH","RH","RH","RH",
"RH","RH","RH","RH")
metric_labels = c("SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"SA","SA","SA","SA",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT",
"CT","CT","CT","CT")
df2plot$region = new_region_labels
df2plot$hemisphere = hemi_labels
df2plot$metric = metric_labels
colors2use = get_ggColorHue(7)
# SA F-stat
d2plot = subset(df2plot, df2plot$metric=="SA" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))
# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot,
display_numbers = round(mat2plot,digits=2),
number_color = "black", fontsize_number = 18,
show_rownames=TRUE,
labels_col = colnames(mat2plot),
color = colorRampPalette(c('light blue','white','red'))(100),
cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
c("#785e00","#c49a00","#ffcc12"),
c("#306800","#53b400","#76ff01"),
c("#007459","#00c094","#0effc8"),
c("#007b9f","#00b6eb","#39d2ff"),
c("#6a3eff","#a58aff","#e0d7ff"),
c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(0.6,1.3), c(0.8,1.15), c(0.75,1.35), c(0.8,1.35),
c(0.8,1.4), c(0.8, 1.15), c(0.7,1.1))
for (ireg in 1:length(sig_reg)){
region = sig_reg[ireg]
reg_name = reg_names[ireg]
metric_name = metric_names[ireg]
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
p1 = p1 + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
print(p1)
tmp_df = d2plot[region,]
if (tmp_df$hemisphere=="LH"){
HEMI = "left"
} else if (tmp_df$hemisphere=="RH"){
HEMI = "right"
}
p = ggseg(.data = tmp_df,
atlas = chenAr,
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere=HEMI) +
scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) +
theme(text = element_text(size=fontSize))
print(p)
}










cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CortexVol"
reg_name = "Total Cortical Volume"
metric_name = "Volume"
ylims = c(450000,700000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p1

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "SAtotal"
reg_name = "Total Surface Area"
metric_name = "SA"
ylims = c(90000,220000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
p1

# SA es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# SA es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# SA es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenAr,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT F-stat
d2plot = subset(df2plot, df2plot$metric=="CT" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))
# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot,
display_numbers = round(mat2plot,digits=2),
number_color = "black", fontsize_number = 18,
show_rownames=TRUE,
labels_col = colnames(mat2plot),
color = colorRampPalette(c('light blue','white','red'))(100),
cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_manual(values = colors2use[1:4]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
# mapping=aes(fill=Fstat),
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_manual(values = colors2use[5:7]) + guides(fill=FALSE) +
# scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
# guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
c("#785e00","#c49a00","#ffcc12"),
c("#306800","#53b400","#76ff01"),
c("#007459","#00c094","#0effc8"),
c("#007b9f","#00b6eb","#39d2ff"),
c("#6a3eff","#a58aff","#e0d7ff"),
c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(-0.4,0.1), c(-0.4,0.4), c(-0.4,0.3), c(0,0.7),
c(-0.3,0.3), c(-0.2,0.5), c(0,0.8))
for (ireg in 1:length(sig_reg)){
region = sig_reg[ireg]
reg_name = reg_names[ireg]
metric_name = metric_names[ireg]
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
p1 = p1 + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) +
xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) +
theme(text = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize),
plot.title = element_text(size=fontSize,hjust=0.5))
print(p1)
tmp_df = d2plot[region,]
if (tmp_df$hemisphere=="LH"){
HEMI = "left"
} else if (tmp_df$hemisphere=="RH"){
HEMI = "right"
}
p = ggseg(.data = tmp_df,
atlas = chenTh,
mapping=aes(fill=label2fill),
position="stacked",
colour="black",
hemisphere=HEMI) +
scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) +
theme(text = element_text(size=fontSize))
print(p)
}




# CT es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.Good_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Good),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf

# CT es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="left") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")
pr = ggseg(.data = df2use,
atlas = chenTh,
mapping=aes(fill=es.TD_vs_Poor),
position="stacked",
colour="black",
hemisphere="right") +
scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
guides(fill = guide_colourbar(nbin = 100)) +
theme(text = element_text(size=fontSize))
pf = (pr)/(pl)
pf
